/*************************************************************
** THROBIT2: Threshold observability bivariate probit model
** Monte Carlo Comparison of Throbit with Boolean Probit
** Skewed to have very few successes
** (c) 2003 by Sanford C. Gordon
**
**
** Department of Politics
** new York University
** Purpose: Determine sample properties of partially observed probit with thresholds.
** output file: mc*.out -- creates Latex-ready tables
*************************************************************/
rndseed 33333;

closeall; 
cls; 

/* Call to the Maxlik library */
library maxlik;
/* Resets global variables to default values */
maxset;

/* Global Parameter Settings */

_max_algorithm = 4;
	                 	@ Optimization Method           @
                       	@ 1 = Steepest Descent          @
                       	@ 2 = BFGS                      @
                       	@ 3 = DFP                       @
                       	@ 4 = Newton-Raphson            @
                       	@ 5 = BHHH                      @
                       	@ 6 = PRCG                      @

_max_CovPar = 1;
						@ Covariance Matrix Method		@
						@ 0 = Not Computed			    @
						@ 1 = Numerical Hessian			@
						@ 2 = Quasi-maximum likelihood  @

_max_LineSearch = 2;
						@ Line Search Method			@
						@ 1 = unit step length			@
						@ 2 = cubic or quadratic		@
						@ 3 = step halving			    @
						@ 4 = Brent's method			@
						@ 5 = BHHH method				@
_max_maxiters = 100;

__output=0;
/*Create symbolic infinity*/
infty=2^10000;


/*******************************************************************************
** Throbit Liklihood Procedure                                                **
*******************************************************************************/

proc throbiti(parms,dta);
    local obs,y,y1,y2,x1,x2,beta1,beta2,tau1,tau2,idx1,idx2,p0,py1,py2,pya,llike;
    /* Data */
    obs = rows(dta);
    y   = dta[.,1];
    y1  = dta[.,2];
    y2  = dta[.,3];
    x1  = ones(obs,1)~dta[.,4:5];
    x2  = ones(obs,1)~dta[.,4]~dta[.,6];
    /* parameters */
    beta1 = parms[1:3];
    beta2 = parms[4:6];
    tau1  = exp(parms[7]);
    tau2  = exp(parms[8]);
    idx1= x1*beta1;
    idx2=x2*beta2;
    p0 = cdfn(-idx1).*cdfn(-idx2); /* This is the probability of a no */
    py1 = cdfn(idx1-tau1).*cdfn(tau2-idx2); /* This is the probability of a yes via mechanism 1 */
    py2 = cdfn(idx2-tau2).*cdfn(tau1-idx1); /* Probability of a yes via mechanism 2 */
    pya = 1-p0-py1-py2; /* Probability of an ambiguous yes */
    /* This is the likelihood function */
    llike = (1-y).*ln(p0)+y.*y1.*ln(py1)+y.*y2.*ln(py2)+y.*(1-y1).*(1-y2).*ln(pya);
    retp(llike);
endp;

/************************************************************************
** Boolean Probit Liklihood Procedure                                  **
************************************************************************/

proc boolprob(parms,dta);
    local y,llike,obs,x1,x2,beta1,beta2,p1,p2,yesprob,noprob;
    obs=rows(dta);
    y=dta[.,1];
    x1  = ones(obs,1)~dta[.,4:5];
    x2  = ones(obs,1)~dta[.,4]~dta[.,6];
    beta1=parms[1:3];
    beta2=parms[4:6];
    p1 = cdfn(x1*beta1);
    p2 = cdfn(x2*beta2);
    yesprob = 1-(1-p1).*(1-p2);  /* Probability of a yes */
    noprob  = 1-yesprob; /* Probability of a no */
    /* This is the likelihood function */
    llike = y.*ln(yesprob)+(1-y).*ln(noprob);
    retp(llike);
endp;

/* Parameter Values Fixed for all simulations */
beta1=1|1|1;
beta2=1|1|1;
tau1 = 0.53;
tau2 = 0.53;
mciter=1000;    /* Number of Monte Carlo simulations */
xoffset = 3.77; /* When drawing Xs, substract this amount to control the outcome balance */


/**************************************************************************************
** First Monte Carlo Routine: Compare throbit and Boolean probit, for sample of 500 **
**************************************************************************************/
/* Generate Xs, which are fixed in repeated samples */
N=500;
xcommon = rndn(N,1);
x12 = rndn(N,1)-xoffset;
x22 = rndn(N,1)-xoffset;
x1 = ones(N,1)~xcommon~x12;
x2 = ones(N,1)~xcommon~x22;

i=1; do until i >mciter;
cls; 
print "N=" n ", i=" i;
/* Generate observed outcomes */
    epsilon1 = rndn(N,1);
    epsilon2 = rndn(N,1);
    zstar1  = x1*beta1+epsilon1;
    zstar2  = x2*beta2+epsilon2;
    y = maxc((zstar1~zstar2)').>0;
    y1= (zstar1.>tau1).*(zstar2.<tau2);
    y2= (zstar2.>tau2).*(zstar1.<tau1);
    let label = y y1 y2 xcommon x12 x22;
    create f0= "c:/gauss50/research/throbit/throbit1" with ^label,0,8;
    call writer(f0,y~y1~y2~xcommon~x12~x22);
    call close(f0);

    dta = "c:/gauss50/research/throbit/throbit1"; /* Tells Maxlik where to look for the data */
    vars = getname(dta);
    

    /* First, run throbit */
    _max_parnames = "beta10"|"beta1c"|"beta12"|"beta20"|"beta2c"|"beta22"|"lntau1"|"lntau2";
    theta0th=0|0|0|0|0|0|0|0;
    {thetath,fth,gth,covth,retth}=maxprt(maxlik(dta,vars,&throbiti,theta0th));

    /* Next, run Boolean probit, using data randomly perturbed around zero */
    _max_parnames = "beta10"|"beta1c"|"beta12"|"beta20"|"beta2c"|"beta22";
    theta0bp=0.01*rndn(6,1);
    {thetabp,fbp,gbp,covbp,retbp}=maxprt(maxlik(dta,vars,&boolprob,theta0bp));

    
    /* Save the results */
    miss6 = {. . . . . .};
    miss8 = {. . . . . . . .};
    if i == 1;
        meany = meanc(y~y1~y2)';
        /* Results for Boolean probit */
        resultbp = thetabp';
        convrgbp = retbp;
        llikebp = fbp*N;
        truthbp = beta1|beta2;
        if retbp == 0;
            msematbp = (thetabp-truthbp)*(thetabp-truthbp)';
        else;
            msematbp = zeros(6,6);
        endif;
        if ismiss(covbp) == 0;
            serrbp = sqrt(diag(covbp))'; 
        else;
            serrbp = miss6;
        endif;
        /* Results for throbit */
        resultth = thetath';
        convrgth = retth;
        lliketh = fth*N;
        truthth = truthbp|ln(tau1)|ln(tau2);
        /* Calculates Mean Squared Error Matrices; #2 is for when Boolean probit fails */
        if retbp == 0;
            msematth = (thetath-truthth)*(thetath-truthth)';
            msemtth2 = zeros(8,8);
        else;
            msemtth2 = (thetath-truthth)*(thetath-truthth)';
            msematth = zeros(8,8);
        endif;
        if ismiss(covth) == 0;
            serrth = sqrt(diag(covth))'; 
        else;
            serrth = miss8;
        endif;
    else;
        meany = meany|meanc(y~y1~y2)';
        resultbp = resultbp|thetabp';
        convrgbp = convrgbp|retbp;
        llikebp = llikebp|fbp*n;
        truthbp = beta1|beta2;
        if retbp ==0; 
            msematbp = msematbp+(thetabp-truthbp)*(thetabp-truthbp)';         
        endif;             
        if ismiss(covbp)==0;
            serrbp=serrbp|sqrt(diag(covbp))';
        else;
            serrbp = serrbp|miss6;
        endif;

        resultth = resultth|thetath';
        convrgth = convrgth|retth;
        lliketh = lliketh|fth*n;
        truthth = truthbp|ln(tau1)|ln(tau2);
        if retbp == 0;
            msematth = msematth+(thetath-truthth)*(thetath-truthth)';
        else;   
            msemtth2 = msemtth2+(thetath-truthth)*(thetath-truthth)';
        endif;
        if ismiss(covth)==0;
            serrth=serrth|sqrt(diag(covth))';
        else;
            serrth = serrth|miss8;
        endif;
    endif;
    i=i+1;
endo;    

sebp500 = meanc(packr(serrbp));
seth500 = meanc(packr(serrth));

keepbp = convrgbp .== 0;  /* This is a vector of zeros and ones; ones correspond to instances where
                             BP had a hard time */
numkeep = sumc(keepbp); /* This is the number of simulations where BP found an answer */
numlose = mciter-numkeep; /* This is the number of simulations where BP failed -- "Hard" Cases */

rmsebp = sqrt(diag(msematbp))/numkeep; /* Calculate root mean squared error for successful cases of BP */
rmseth = sqrt(diag(msematth))/numkeep; /* Calculate root mean squared error for throbit estimates 
                                          corresponding to cases in which BP succeeded */
if numlose == 0;
    rmseth2 = zeros(8,1);
else;
    rmseth2= sqrt(diag(msemtth2))/numlose; /* Calculate rmse for throbit when BP failed */
endif;

resbpk = delif(resultbp,1-keepbp); /*Results for successful Boolean probit */
resthk = delif(resultth,1-keepbp); /*Results for throbit when BP was successful */
resthl = delif(resultth,keepbp);   /*Results for throbit when BP failed */

/* Print output table */ 
amper4= "&"$|"&"$|"&"$|"&";
slashes4  = "\\\\"$|"\\\\"$|"\\\\"$|"\\\\";
caption = "Truth"|"Mean Boolean"|"Mean Throbit"|"Rel. Eff.";
meanbp = meanc(resbpk);
meanth = meanc(resthk);
releff = rmsebp./(rmseth[1:6]);
b10 = truthbp[1]|meanbp[1]|meanth[1]|releff[1];
b1c = truthbp[2]|meanbp[2]|meanth[2]|releff[2];
b12 = truthbp[3]|meanbp[3]|meanth[3]|releff[3];
b20 = truthbp[4]|meanbp[4]|meanth[4]|releff[4]; 
b2c = truthbp[5]|meanbp[5]|meanth[5]|releff[5];
b22 = truthbp[6]|meanbp[6]|meanth[6]|releff[6];
lnt1 = meanth[7];
lnt2 = meanth[8];
lnt1t = ln(tau1);
lnt2t = ln(tau2);

b10s = "" $+ftocv(b10,0,2);
b1cs = "" $+ftocv(b1c,0,2);
b12s = "" $+ftocv(b12,0,2);
b20s = "" $+ftocv(b20,0,2);
b2cs = "" $+ftocv(b2c,0,2);
b22s = "" $+ftocv(b22,0,2);
lnt1s = "" $+ftocv(lnt1t,0,2)|"--"|"" $+ftocv(lnt1,0,2)|"--";
lnt2s = "" $+ftocv(lnt2t,0,2)|"--"|"" $+ftocv(lnt2,0,2)|"--";

resmat = caption$~amper4$~b10s$~amper4$~b1cs$~amper4$~b12s$~amper4$~b20s$~amper4$~b2cs$~amper4$~b22s$~amper4$~lnt1s$~amper4$~lnt2s$~slashes4;
output file = "c:/gauss50/research/throbit/mc2.out" reset;
print"\\begin{table}";
print"\\caption{\\label{mc2} Monte Carlo Comparison of Boolean Probit and Throbit, Skewed Outcomes}";
print"\\begin{center}";
print"\\begin{tabular}{lcccccccc}";
print"\\hline \\hline";
print"&$\\beta_{10}$&$\\beta_{1c}$&$\\beta_{12}$&$\\beta_{20}$&$\\beta_{2c}$&$\\beta_{22}$&$\\ln(\\tau_1)$&$\\ln(\\tau_2)$ \\\\ ";

print "\\hline";
print "\\multicolumn{1}{l}{} &  \\multicolumn{8}{c}{\\emph{N=500}} \\\\ \\hline";
resmat;
print "\\hline";
print "\\multicolumn{9}{l}{Kept simulations: " ""$+ftocv(numkeep,0,0)$+"}\\\\";
output off;
/* Print Hard Cases Table */

amper3= "&"$|"&"$|"&";
slashes3  = "\\\\"$|"\\\\"$|"\\\\";
caption = "Truth"|"Mean Throbit"|"RMSE";
if ismiss(resthl)==1;
    meanthl = zeros(8,1);
else;
    meanthl = meanc(resthl);
endif;
b10 = truthth[1]|meanthl[1]|rmseth2[1];
b1c = truthth[2]|meanthl[2]|rmseth2[2];
b12 = truthth[3]|meanthl[3]|rmseth2[3];
b20 = truthth[4]|meanthl[4]|rmseth2[4];
b2c = truthth[5]|meanthl[5]|rmseth2[5];
b22 = truthth[6]|meanthl[6]|rmseth2[6];
lnt1 = truthth[7]|meanthl[7]|rmseth2[7];
lnt2 = truthth[8]|meanthl[8]|rmseth2[8];

b10s = "" $+ftocv(b10,0,2);
b1cs = "" $+ftocv(b1c,0,2);
b12s = "" $+ftocv(b12,0,2);
b20s = "" $+ftocv(b20,0,2);
b2cs = "" $+ftocv(b2c,0,2);
b22s = "" $+ftocv(b22,0,2);
lnt1s = "" $+ftocv(lnt1,0,2);
lnt2s = "" $+ftocv(lnt2,0,2);

resmat2 = caption$~amper3$~b10s$~amper3$~b1cs$~amper3$~b12s$~amper3$~b20s$~amper3$~b2cs$~amper3$~b22s$~amper3$~lnt1s$~amper3$~lnt2s$~slashes3;

output file = "c:/gauss50/research/throbit/hardcase.out" reset;
print"\\begin{table}";
print"\\caption{\\label{hardcase} Monte Carlo Simulation of Throbit, ``Hard'' Cases}";
print"\\begin{center}";
print"\\begin{tabular}{lcccccccc}";
print"\\hline \\hline";
print"&$\\beta_{10}$&$\\beta_{1c}$&$\\beta_{12}$&$\\beta_{20}$&$\\beta_{2c}$&$\\beta_{22}$&$\\ln(\\tau_1)$&$\\ln(\\tau_2)$ \\\\ ";
print "\\hline";
print "\\multicolumn{1}{l}{} &  \\multicolumn{8}{c}{\\emph{N=500}} \\\\ \\hline";
resmat2;
print "\\hline";
print "\\multicolumn{9}{l}{Number of Hard Cases: " ""$+ftocv(numlose,0,0)$+"}\\\\";
output off;

output file = "c:/gauss50/research/throbit/diag2.out" on;
print "I have finished with N=500, moving on to N=1000";
output off;



/**************************************************************************************
** Second Monte Carlo Routine: Compare throbit and Boolean probit, for sample of 1000 **
**************************************************************************************/
/* See above for more detailed comments -- everything is the same except different 
   sample size and table output doesn't have frontmatter */

/* Generate Xs, which are fixed in repeated samples */
N=1000;
xcommon = rndn(N,1);
x12 = rndn(N,1)-xoffset;
x22 = rndn(N,1)-xoffset;
x1 = ones(N,1)~xcommon~x12;
x2 = ones(N,1)~xcommon~x22;

i=1; do until i >mciter;
cls; 
print "N=" n ", i=" i;
/* Generate observed outcomes */
    epsilon1 = rndn(N,1);
    epsilon2 = rndn(N,1);
    zstar1  = x1*beta1+epsilon1;
    zstar2  = x2*beta2+epsilon2;
    y = maxc((zstar1~zstar2)').>0;
    y1= (zstar1.>tau1).*(zstar2.<tau2);
    y2= (zstar2.>tau2).*(zstar1.<tau1);
    let label = y y1 y2 xcommon x12 x22;
    create f0= "c:/gauss50/research/throbit/throbit1" with ^label,0,8;
    call writer(f0,y~y1~y2~xcommon~x12~x22);
    call close(f0);

    dta = "c:/gauss50/research/throbit/throbit1"; /* Tells Maxlik where to look for the data */
    vars = getname(dta);
    

    /* First, run throbit */
    _max_parnames = "beta10"|"beta1c"|"beta12"|"beta20"|"beta2c"|"beta22"|"lntau1"|"lntau2";
    theta0th=0|0|0|0|0|0|0|0;
    {thetath,fth,gth,covth,retth}=maxprt(maxlik(dta,vars,&throbiti,theta0th));

    /* Next, run Boolean probit, using throbit parameters as starts */
    _max_parnames = "beta10"|"beta1c"|"beta12"|"beta20"|"beta2c"|"beta22";
    theta0bp=0.01*rndn(6,1);
    {thetabp,fbp,gbp,covbp,retbp}=maxprt(maxlik(dta,vars,&boolprob,theta0bp));
    
    /* Save the results */
    miss6 = {. . . . . .};
    miss8 = {. . . . . . . .};
    if i == 1;
        meany = meanc(y~y1~y2)';
        /* Results for Boolean probit */
        resultbp = thetabp';
        convrgbp = retbp;
        llikebp = fbp*N;
        truthbp = beta1|beta2;
        if retbp == 0;
            msematbp = (thetabp-truthbp)*(thetabp-truthbp)';
        else;
            msematbp = zeros(6,6);
        endif;
        if ismiss(covbp) == 0;
            serrbp = sqrt(diag(covbp))'; 
        else;
            serrbp = miss6;
        endif;
        /* Results for throbit */
        resultth = thetath';
        convrgth = retth;
        lliketh = fth*N;
        truthth = truthbp|ln(tau1)|ln(tau2);
        if retbp == 0;
            msematth = (thetath-truthth)*(thetath-truthth)';
            msemtth2 = zeros(8,8);
        else;
            msemtth2 = (thetath-truthth)*(thetath-truthth)';
            msematth = zeros(8,8);
        endif;
        if ismiss(covth) == 0;
            serrth = sqrt(diag(covth))'; 
        else;
            serrth = miss8;
        endif;
    else;
        meany = meany|meanc(y~y1~y2)';
        resultbp = resultbp|thetabp';
        convrgbp = convrgbp|retbp;
        llikebp = llikebp|fbp*n;
        truthbp = beta1|beta2;
        if retbp ==0; 
            msematbp = msematbp+(thetabp-truthbp)*(thetabp-truthbp)';
        endif;             
        if ismiss(covbp)==0;
            serrbp=serrbp|sqrt(diag(covbp))';
        else;
            serrbp = serrbp|miss6;
        endif;

        resultth = resultth|thetath';
        convrgth = convrgth|retth;
        lliketh = lliketh|fth*n;
        truthth = truthbp|ln(tau1)|ln(tau2);
        if retbp == 0;
            msematth = msematth+(thetath-truthth)*(thetath-truthth)';
        else;   
            msemtth2 = msemtth2+(thetath-truthth)*(thetath-truthth)';
        endif;
        if ismiss(covth)==0;
            serrth=serrth|sqrt(diag(covth))';
        else;
            serrth = serrth|miss8;
        endif;
    endif;
    i=i+1;
endo;    

sebp500 = meanc(packr(serrbp));
seth500 = meanc(packr(serrth));

keepbp = convrgbp .== 0;
numkeep = sumc(keepbp);
numlose = mciter-numkeep;

rmsebp = sqrt(diag(msematbp))/numkeep;
rmseth = sqrt(diag(msematth))/numkeep;
if numlose == 0;
    rmseth2 = zeros(8,1);
else;
    rmseth2= sqrt(diag(msemtth2))/numlose;
endif;

resbpk = delif(resultbp,1-keepbp);
resthk = delif(resultth,1-keepbp);
resthl = delif(resultth,keepbp);

/* Print table */ 
amper4= "&"$|"&"$|"&"$|"&";
slashes4  = "\\\\"$|"\\\\"$|"\\\\"$|"\\\\";
caption = "Truth"|"Mean Boolean"|"Mean Throbit"|"Rel. Eff.";
meanbp = meanc(resbpk);
meanth = meanc(resthk);
releff = rmsebp./(rmseth[1:6]);
b10 = truthbp[1]|meanbp[1]|meanth[1]|releff[1];
b1c = truthbp[2]|meanbp[2]|meanth[2]|releff[2];
b12 = truthbp[3]|meanbp[3]|meanth[3]|releff[3];
b20 = truthbp[4]|meanbp[4]|meanth[4]|releff[4]; 
b2c = truthbp[5]|meanbp[5]|meanth[5]|releff[5];
b22 = truthbp[6]|meanbp[6]|meanth[6]|releff[6];
lnt1 = meanth[7];
lnt2 = meanth[8];
lnt1t = ln(tau1);
lnt2t = ln(tau2);

b10s = "" $+ftocv(b10,0,2);
b1cs = "" $+ftocv(b1c,0,2);
b12s = "" $+ftocv(b12,0,2);
b20s = "" $+ftocv(b20,0,2);
b2cs = "" $+ftocv(b2c,0,2);
b22s = "" $+ftocv(b22,0,2);
lnt1s = "" $+ftocv(lnt1t,0,2)|"--"|"" $+ftocv(lnt1,0,2)|"--";
lnt2s = "" $+ftocv(lnt2t,0,2)|"--"|"" $+ftocv(lnt2,0,2)|"--";

resmat = caption$~amper4$~b10s$~amper4$~b1cs$~amper4$~b12s$~amper4$~b20s$~amper4$~b2cs$~amper4$~b22s$~amper4$~lnt1s$~amper4$~lnt2s$~slashes4;
output file = "c:/gauss50/research/throbit/mc2.out" on;
print "\\hline";
print "\\multicolumn{1}{l}{} &  \\multicolumn{8}{c}{\\emph{N=1,000}} \\\\ \\hline";
resmat;
print "\\hline";
print "\\multicolumn{9}{l}{Kept simulations: " ""$+ftocv(numkeep,0,0)$+"}\\\\";
output off;
/* Print Hard Cases Table */

amper3= "&"$|"&"$|"&";
slashes3  = "\\\\"$|"\\\\"$|"\\\\";
caption = "Truth"|"Mean Throbit"|"RMSE";
if ismiss(resthl)==1;
    meanthl = zeros(8,1);
else;
    meanthl = meanc(resthl);
endif;
b10 = truthth[1]|meanthl[1]|rmseth2[1];
b1c = truthth[2]|meanthl[2]|rmseth2[2];
b12 = truthth[3]|meanthl[3]|rmseth2[3];
b20 = truthth[4]|meanthl[4]|rmseth2[4];
b2c = truthth[5]|meanthl[5]|rmseth2[5];
b22 = truthth[6]|meanthl[6]|rmseth2[6];
lnt1 = truthth[7]|meanthl[7]|rmseth2[7];
lnt2 = truthth[8]|meanthl[8]|rmseth2[8];

b10s = "" $+ftocv(b10,0,2);
b1cs = "" $+ftocv(b1c,0,2);
b12s = "" $+ftocv(b12,0,2);
b20s = "" $+ftocv(b20,0,2);
b2cs = "" $+ftocv(b2c,0,2);
b22s = "" $+ftocv(b22,0,2);
lnt1s = "" $+ftocv(lnt1,0,2);
lnt2s = "" $+ftocv(lnt2,0,2);

resmat2 = caption$~amper3$~b10s$~amper3$~b1cs$~amper3$~b12s$~amper3$~b20s$~amper3$~b2cs$~amper3$~b22s$~amper3$~lnt1s$~amper3$~lnt2s$~slashes3;

output file = "c:/gauss50/research/throbit/hardcase.out" on;
print "\\hline";
print "\\multicolumn{1}{l}{} &  \\multicolumn{8}{c}{\\emph{N=1,000}} \\\\ \\hline";
resmat2;
print "\\hline";
print "\\multicolumn{9}{l}{Number of Hard Cases: " ""$+ftocv(numlose,0,0)$+"}\\\\";
output off;

output file = "c:/gauss50/research/throbit/diag2.out" on;
print "I have finished with N=1000, moving on to N=5000";
output off;


/**************************************************************************************
** Third Monte Carlo Routine: Compare throbit and Boolean probit, for sample of 5000 **
**************************************************************************************/
/* See above for more detailed comments -- everything is the same except different 
   sample size and table output doesn't have frontmatter */

/* Generate Xs, which are fixed in repeated samples */
N=5000;
xcommon = rndn(N,1);
x12 = rndn(N,1)-xoffset;
x22 = rndn(N,1)-xoffset;
x1 = ones(N,1)~xcommon~x12;
x2 = ones(N,1)~xcommon~x22;

i=1; do until i >mciter;
cls; 
print "N=" n ", i=" i;
/* Generate observed outcomes */
    epsilon1 = rndn(N,1);
    epsilon2 = rndn(N,1);
    zstar1  = x1*beta1+epsilon1;
    zstar2  = x2*beta2+epsilon2;
    y = maxc((zstar1~zstar2)').>0;
    y1= (zstar1.>tau1).*(zstar2.<tau2);
    y2= (zstar2.>tau2).*(zstar1.<tau1);
    let label = y y1 y2 xcommon x12 x22;
    create f0= "c:/gauss50/research/throbit/throbit1" with ^label,0,8;
    call writer(f0,y~y1~y2~xcommon~x12~x22);
    call close(f0);

    dta = "c:/gauss50/research/throbit/throbit1"; /* Tells Maxlik where to look for the data */
    vars = getname(dta);
    

    /* First, run throbit */
    _max_parnames = "beta10"|"beta1c"|"beta12"|"beta20"|"beta2c"|"beta22"|"lntau1"|"lntau2";
    theta0th=0|0|0|0|0|0|0|0;
    {thetath,fth,gth,covth,retth}=maxprt(maxlik(dta,vars,&throbiti,theta0th));

    /* Next, run Boolean probit, using throbit parameters as starts */
    _max_parnames = "beta10"|"beta1c"|"beta12"|"beta20"|"beta2c"|"beta22";
    theta0bp=0.01*rndn(6,1);
    {thetabp,fbp,gbp,covbp,retbp}=maxprt(maxlik(dta,vars,&boolprob,theta0bp));

    
    /* Save the results */
    miss6 = {. . . . . .};
    miss8 = {. . . . . . . .};
    if i == 1;
        meany = meanc(y~y1~y2)';
        /* Results for Boolean probit */
        resultbp = thetabp';
        convrgbp = retbp;
        llikebp = fbp*N;
        truthbp = beta1|beta2;
        if retbp == 0;
            msematbp = (thetabp-truthbp)*(thetabp-truthbp)';
        else;
            msematbp = zeros(6,6);
        endif;
        if ismiss(covbp) == 0;
            serrbp = sqrt(diag(covbp))'; 
        else;
            serrbp = miss6;
        endif;
        /* Results for throbit */
        resultth = thetath';
        convrgth = retth;
        lliketh = fth*N;
        truthth = truthbp|ln(tau1)|ln(tau2);
        if retbp == 0;
            msematth = (thetath-truthth)*(thetath-truthth)';
            msemtth2 = zeros(8,8);
        else;
            msemtth2 = (thetath-truthth)*(thetath-truthth)';
            msematth = zeros(8,8);
        endif;
        if ismiss(covth) == 0;
            serrth = sqrt(diag(covth))'; 
        else;
            serrth = miss8;
        endif;
    else;
        meany = meany|meanc(y~y1~y2)';
        resultbp = resultbp|thetabp';
        convrgbp = convrgbp|retbp;
        llikebp = llikebp|fbp*n;
        truthbp = beta1|beta2;
        if retbp ==0; 
            msematbp = msematbp+(thetabp-truthbp)*(thetabp-truthbp)';
        endif;             
        if ismiss(covbp)==0;
            serrbp=serrbp|sqrt(diag(covbp))';
        else;
            serrbp = serrbp|miss6;
        endif;

        resultth = resultth|thetath';
        convrgth = convrgth|retth;
        lliketh = lliketh|fth*n;
        truthth = truthbp|ln(tau1)|ln(tau2);
        if retbp == 0;
            msematth = msematth+(thetath-truthth)*(thetath-truthth)';
        else;   
            msemtth2 = msemtth2+(thetath-truthth)*(thetath-truthth)';
        endif;
        if ismiss(covth)==0;
            serrth=serrth|sqrt(diag(covth))';
        else;
            serrth = serrth|miss8;
        endif;
    endif;
    i=i+1;
endo;    

sebp500 = meanc(packr(serrbp));
seth500 = meanc(packr(serrth));

keepbp = convrgbp .== 0;
numkeep = sumc(keepbp);
numlose = mciter-numkeep;

rmsebp = sqrt(diag(msematbp))/numkeep;
rmseth = sqrt(diag(msematth))/numkeep;
if numlose == 0;
    rmseth2 = zeros(8,1);
else;
    rmseth2= sqrt(diag(msemtth2))/numlose;
endif;

resbpk = delif(resultbp,1-keepbp);
resthk = delif(resultth,1-keepbp);
resthl = delif(resultth,keepbp);

/* Print table */ 
amper4= "&"$|"&"$|"&"$|"&";
slashes4  = "\\\\"$|"\\\\"$|"\\\\"$|"\\\\";
caption = "Truth"|"Mean Boolean"|"Mean Throbit"|"Rel. Eff.";
meanbp = meanc(resbpk);
meanth = meanc(resthk);
releff = rmsebp./(rmseth[1:6]);
b10 = truthbp[1]|meanbp[1]|meanth[1]|releff[1];
b1c = truthbp[2]|meanbp[2]|meanth[2]|releff[2];
b12 = truthbp[3]|meanbp[3]|meanth[3]|releff[3];
b20 = truthbp[4]|meanbp[4]|meanth[4]|releff[4]; 
b2c = truthbp[5]|meanbp[5]|meanth[5]|releff[5];
b22 = truthbp[6]|meanbp[6]|meanth[6]|releff[6];
lnt1 = meanth[7];
lnt2 = meanth[8];
lnt1t = ln(tau1);
lnt2t = ln(tau2);

b10s = "" $+ftocv(b10,0,2);
b1cs = "" $+ftocv(b1c,0,2);
b12s = "" $+ftocv(b12,0,2);
b20s = "" $+ftocv(b20,0,2);
b2cs = "" $+ftocv(b2c,0,2);
b22s = "" $+ftocv(b22,0,2);
lnt1s = "" $+ftocv(lnt1t,0,2)|"--"|"" $+ftocv(lnt1,0,2)|"--";
lnt2s = "" $+ftocv(lnt2t,0,2)|"--"|"" $+ftocv(lnt2,0,2)|"--";

resmat = caption$~amper4$~b10s$~amper4$~b1cs$~amper4$~b12s$~amper4$~b20s$~amper4$~b2cs$~amper4$~b22s$~amper4$~lnt1s$~amper4$~lnt2s$~slashes4;
output file = "c:/gauss50/research/throbit/mc2.out" on;
print "\\hline";
print "\\multicolumn{1}{l}{} &  \\multicolumn{8}{c}{\\emph{N=1,000}} \\\\ \\hline";
resmat;
print "\\hline";
print "\\multicolumn{9}{l}{Kept simulations: " ""$+ftocv(numkeep,0,0)$+"}\\\\";
print " \\hline \\hline ";
print "\\end{tabular}";
print "\\end{center}";
print "\\end{table}";
output off;
/* Print Hard Cases Table */

amper3= "&"$|"&"$|"&";
slashes3  = "\\\\"$|"\\\\"$|"\\\\";
caption = "Truth"|"Mean Throbit"|"RMSE";
if ismiss(resthl)==1;
    meanthl = zeros(8,1);
else;
    meanthl = meanc(resthl);
endif;

b10 = truthth[1]|meanthl[1]|rmseth2[1];
b1c = truthth[2]|meanthl[2]|rmseth2[2];
b12 = truthth[3]|meanthl[3]|rmseth2[3];
b20 = truthth[4]|meanthl[4]|rmseth2[4];
b2c = truthth[5]|meanthl[5]|rmseth2[5];
b22 = truthth[6]|meanthl[6]|rmseth2[6];
lnt1 = truthth[7]|meanthl[7]|rmseth2[7];
lnt2 = truthth[8]|meanthl[8]|rmseth2[8];

b10s = "" $+ftocv(b10,0,2);
b1cs = "" $+ftocv(b1c,0,2);
b12s = "" $+ftocv(b12,0,2);
b20s = "" $+ftocv(b20,0,2);
b2cs = "" $+ftocv(b2c,0,2);
b22s = "" $+ftocv(b22,0,2);
lnt1s = "" $+ftocv(lnt1,0,2);
lnt2s = "" $+ftocv(lnt2,0,2);

resmat2 = caption$~amper3$~b10s$~amper3$~b1cs$~amper3$~b12s$~amper3$~b20s$~amper3$~b2cs$~amper3$~b22s$~amper3$~lnt1s$~amper3$~lnt2s$~slashes3;

output file = "c:/gauss50/research/throbit/hardcase.out" on;
print "\\hline";
print "\\multicolumn{1}{l}{} &  \\multicolumn{8}{c}{\\emph{N=1,000}} \\\\ \\hline";
resmat2;
print "\\hline";
print "\\multicolumn{9}{l}{Number of Hard Cases: " ""$+ftocv(numlose,0,0)$+"}\\\\";
print " \\hline \\hline ";
print "\\end{tabular}";
print "\\end{center}";
print "\\end{table}";






output off;

output file = "c:/gauss50/research/throbit/diag2.out" on;
print "I have finished with N=5000.";
output off;



